Crystallization of amorphous silicon induced by mechanical shear deformations 
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We have investigated the response of amorphous sihcon (a-Si), in particular crystaUization, to 
external mechanical shear deformations using classical molecular dynamics (MD) simulations and the 
empirical Environment Dependent Inter-atomic Potential (EDIP) [Phys. Rev. B 56, 8542 (1997)]. 
In agreement with previous results we find that, at low shear velocity and low temperature, shear 
deformations increase disorder and defect density. At high temperatures, however, the deformations 
are found to induce crystallization, demonstrating a dynamical transition associated with both shear 
rate and temperature. The properties of a- Si under shear deformations and the extent at which 
the system crystallizes are analyzed in terms of the potential energy difference (PED) between the 
sheared and non-sheared material, as well as the fraction of defects and the number of particles that 
possess a crystalline environment. 
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I. INTRODUCTION 

Amorphous silicon (a-Si) is widely used in the elec- 
tronic and semiconductor industries. From a theoretical 
point of view, because of its simplicity, it has become the 
archetype of covalent amorphous systems that include sil- 
ica and chalcogenide glasses. While these materials have 
been extensively studied for almost four decades, their 
behavior under external forces has received significant 
attention only recently yet, many questions remain 
unanswered. For example, while the plasticity of a-Si is 
attributed to the presence of liquid-like particles associ- 
ated with 5-fold coordinated atoms, ^ ^ a complete picture 
of the deformation effects on disordered and amorphous 
materials is far from being complete. It has been shown 
that, under external forces, amorphous materials may 
crystallize,!^^"^ providing insight into the nucleation pro- 
cess for these covalent materials which is only partially 
understood. Indeed, crystallization is often studied 
using amorphous-crystal or liquid-cryst al interf acej^^^^ 
and this is in particular the case for SifHEmnHUl Crys- 
tallization may also result from the application of exter- 
nal forces such as mechanical shear deformation^^^M^ or 
magnetic fields.'^ 

The properties of bulk crystalline silicon (c-Si) have 
been investigated by means of molecular dynamics (MD) 
simulations using both classical potential^^ and ah- 
initio methods It has been demonstrated that, un- 
der the effect of shear, the interatomic bonds lose their 
covalent character until a metallic state is established 
in the shear directionpl upon inverting the deforma- 
tion, however, the perfect diamond structure is recov- 
ered. In the case of disordered materials, it has been 
shown that small strains bring these systems to deeper 
energy minima in the glassy state. Under deformations, 
these materials might escape the minima through the 
high energy barrier and visit other minima of the energy 
landscape giving rise to new local rearrangements of the 



amorphous structure. These features depend on strain, 
shear rate and t emperature at which the deformations 
are applied. I^^^*^ 

Strain can also play a role in re-ordering the network. 
Lee et a/.,'^for example, have investigated deformation 
and grain growth in partially crystallized nickel by means 
of MD simulations and showed that shear deformations 
can enhance crystallization in am orphous materials. Re- 
cently, Mokshin and Barrat^^ have demonstrated that 
shearing an initially amorphous system leads to an in- 
crease of the nano-crystalline ordering. They examined 
both a one-component Lennard- Jones system and the 
one-component short-range Dzugutov potential. As the 
disordered phase of these single-component close-packed 
systems is known to be very unstable (even though the 
Dzugutov model does better in this respect^^), however, 
it is not clear whether this shear-induced crystallization 
is generic. Talati et al.J using both the Tersofi"^^ po- 
tential and various versions of the Stillinger-Webei^^ po- 
tential for silicon, found no evidence of ordering from 
the amorphous state in the range of shear velocities 
and temperatures they considered. The crystallization 
of computer-simulated silicon is, however, known to be 
potential dependent, and both standard Tersoff and 
the Stillinger- Weber potentials tend to over-stabilize the 
liquid-like environments in the amorphous phase. Given 
the importance of shear-induced crystallization, further 
study of this covalent system is clearly warranted. 

In view of this situation, we have carried out a detailed 
MD study of the properties of a-Si under low-velocity 
shear deformations, describing the energetics of the sili- 
con atoms in terms of the Envi ronment Dependent Inter- 
atomic Potential (EDIP). The present study follows 
up on a previous investigation of the effects of shear ve- 
locity and temperature on amorphous silicon. In partic- 
ular, it deals with the extent of the distribution of defects 
when a low shear velocity is imposed on the system or 
when the deformations are applied at high temperature. 
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Our previous study was concerned with the effects of 
shear velocity and temperature on amorphous sihcon; it 
was found that the impact of an externahy apphed strain 
can be almost fully compensated by increasing the tem- 
perature, allowing the system to respond more rapidly 
to the deformation. In this article we show that that, 
depending on shear velocity and temperature, a- Si can 
either crystallize or remain in its glassy state, following a 
more complex kinetic path than the Lennard- Jones sys- 
tem but consistent with the "universality" proposed by 
Mokshin and Barrat.l^^ 



II. COMPUTER MODEL 

As mentioned above, our calculations were performed 
within t he fra mework of MD and the EDIP classical 
potential l^^l^ This potential, which was fitted to vari- 
ous zero-temperature bulk phases and defect structures, 
reproduces accurately the structure and the dynamics 
of a- Si as well as different crystalline structures. It 
possesses remarkable transferability for zero-temperature 
properties, including elas tic co nstants, bulk crystal struc- 
tures and point defect s.*^^^ This potential has been 
used with succ ess t o study ion-beam induced amor- 
phization of C-SPEI 

as well as crystalliza tion o f a- Si at 
the amorphous-crystalline silicon interface.'^^'^ All sim- 
ulations were performed using the massively parallel 
MD package LAMMPS developed by Sandia National 
Laboratories^ 

The method adopted to apply the mechanical shear de- 
formations has been discussed in detail in Ref.^^ Three 
regions are defined, as illustrated in Fig. [ij an upper 
and a lower walls (perpendicular to the ?/-axis) used to 
apply the mechanical shear deformations, and a central 
region, with unconstrained mobile particles. The parti- 
cles in the lower wall are fixed while those in the upper 
wall move as a whole in the shear direction (x) at fixed 
velocity Vg (shear velocity). The distance between the 
two walls (c:^ 42 A) is kept constant during the simula- 
tions. The equations of motion for the mobile particles 
are integrated using the velocity- Verlet algorithm with a 
time step of 1 fs. Periodic boundary conditions are im- 
posed in the x and z directions, while they are fixed in 
the y direction. The thickness of the walls {c::^^ A) is cho- 
sen such as they each contain 1000 particles; the mobile 
region thus contains 6000 particles. The walls are made 
out of material from the a- Si matrix, so that their struc- 
ture is similar to that of the bulk before shear is applied; 
this choice ensures that no bias towards crystal growth 
or structural modifications can be induced by the wall 
layers. 

To follow the system's response to shear deformations, 
we compute the radial distribution function (RDF), the 
potential energy difference (PED) and the fraction of de- 
fects (5-fold coordinated atoms) between the deformed 
system and the initial non-sheared model.!^ To follow a 
possible crystallization process, we employ a topological 




FIG. 1: (color online) A typical snapshot of the (8000- atom) 
a-Si model; the mobile particles are in the center region be- 
tween two parallel walls (see text for details). 



description based on the ring structure of the two small 
crystalline building blocks that characterize the diamond 
and the hexagonal structures.*^ 



III. SHEAR DEFORMATION RESULTS 
A. Effects of shearing as a function of temperature 

We have shown, in a previous study,^^ that high shear 
rates (relative to the temperature-dependent relaxation 
rate) induce strong and inhomogeneous local strains in 
the plastic regime; to allow the system to adapt to the 
shear deformations, it is therefore necessary to drive it as 
slowly as possible. On the basis of this work, we select a 
shear velocity of 10~^ A/ps, the slowest rate that allows 
the system to reach the plastic regime at all temperatures 
studied in the time scale accessible by MD simulations. 

Simulations were performed at four different tempera- 
tures (300, 600, 900 and 1000 K) and were stopped when 
the total strain reached 12 %, which is sufficient to in- 
duce plastic deformations of amorphous silicon. The 
total simulation time of each run is about 500 ns. Fig- 
ure m shows the evolution of the PED between sheared 
and non-sheared system as a function of imposed strain 
for the different temperatures considered. While a sys- 
tematic increase of the PED was observed at all temper- 
atures when the deformations are driven with high shear 
velocity,^^ Fig. [2] presents a more complex picture with 
three well-defined regimes: an overall energy increase for 
the two lowest temperatures, slight relaxation at 900 K, 
and clear energy decrease at 1000 K. 

More precisely, for the two lowest temperatures, 300 
and 600 K, the potential energy increases by 0.02 and 
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FIG. 2: (color online) PED between sheared and non-sheared 
systems at temperatures 300, 600, 900, and 1000 K. The arrow 
indicates increasing temperature. 



0.004 eV/atom, respectively, from the initial equilibrium 
state to the steady plastic regime. This behavior is 
characteristic of strain-induced disordering, in agreement 
with the results discussed by us in Ref. |33 At 900 K, 
there is no peak associated with the elastic to plastic 
transition; we observe rather a drop in energy at low 
strain to a plateau 0.005 eV/atom below the initial value. 
Shearing thus allows the system to reach a lower-energy 
state — still amorphous — as had been observed, for 
example, by Isner and Lacks. At 1000 K, finally, the 
system's response to the deformation is qualitatively dif- 
ferent: the energy falls almost linearly up to a strain of 
about 10 %, reaching a value about 0.06 eV/atom below 
the initial state. The energy continues to fall, albeit at 
a slower rate, for strains up to 12 %. Clearly, the en- 
ergy change resulting from the application of shear is the 
manifestation of strong structural modifications in the 
network. 

To understand the microscopic origin of these changes, 
it is useful to examine the behavior of some structural 
parameters. We present, in Fig. [3) the variation as a 
function of strain of the proportion of 4- and 5-fold coor- 
dinated atoms for the four temperatures considered. We 
have demonstrated, in a previous article, that the in- 
crease in potential energy of a-Si resulting from shear is 
related to the increase in the number of defects, mainly 
5-fold coordinated atoms that are considered as "liquid- 
like". As suggested in Refs. 5 and 6, these defective 
liquid-like atoms are directly associated with plastic de- 
formations in amorphous silicon. 

Fig. |3] demonstrates, here also, a direct correlation be- 
tween changes in the potential energy and defects cre- 
ated by the shear deformation. Our simulations show 
that the number of 5-fold atoms (liquid-like] atoms first 
increases as a function of shear at low temperatures, but 
decreases at 900 K and 1000 K. At equilibrium (zero 
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FIG. 3: (color online) 4- (top panel) and 5-fold (bottom panel) 
coordinated atoms as a function of strain at temperatures 300, 
600, 900, and 1000 K. A cutoff of 2.8 A is used for d efining 
which atom are included in the nearest-neighbour count. The 
arrow indicates increasing temperature. 



strain), the fraction of 5-fold defects increases with tem- 
perature, from 1.5 % at 300 K to 4.9 % at 1000 K, as can 
be seen in Fig.[3| Upon shearing, in contrast, the system 
crosses over to the opposite situation: the relative num- 
ber of defects decreases as temperature increases. Thus, 
from almost 7 % at 300 K for maximum strain — the 
system is unable to relax the defects resulting from the 
shear — , it drops to 4 %, at 600K and to 1.6 % at 1000 K. 
As with the PED, the 900 K simulation shows the least 
impact under shear, with the proportion of defects drop- 
ping from 3.5 % to about about 2.5%atal2% strain. 

To characterize the response of amorphous silicon to 
shear deformations, we show in Fig. [4] a 3D representa- 
tion of the distribution of defects as a function of strain 
and layer index (distance from the edges of the simula- 
tion cell — see caption). This figure reveals that, in the 
plastic regime, the defects are mostly concentrated in a 
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few layers at 300 K, but that this distribution widens 
as temperature is raised and becomes uniform at 900 K. 
While the defect concentration is similar between 900 
and 1000 K, we note that the distribution are qualita- 
tively different and defects are concentrated at the in- 
terface with the rigid walls for the highest temperature 
(Fig.[4](d)). The decrease of the defect fraction at 900 K 
is consistent with the small relaxation observed in the 
PED. At 1000 K, the small proportion and flat distri- 
bution of defects over several layers indicate significant 
changes in the structure of the sheared material. As we 
discuss below, this behavior is associated with crystalliza- 
tion spreading over a large region about the center of the 
system. Since the walls are frozen in the a-Si structure, 
the extent of crystallization is limited by the presence of 
the walls and liquid-like defects are likely to occur at the 
amorphous-crystal interface. 

To obtain a first characterization of the structural sig- 
nificance of these changes, we examine the radial distri- 
bution function (RDF) of these four models at maximum 
strain (12 %) and compare it with that of unstrained a- 
Si at 300 K (Fig. [5|. At 12 % strain, the energy and 
the defect density are saturated in all cases and are thus 
representative of the steady state plastic regime. 

For temperatures 300, 600, and 900 K, the RDFs are 
very similar and resemble closely that of relaxed a-Si. 
The positions of the first- and second-neighbor peaks co- 
incide, as evident in the inset. The main manifestation of 
shearing is the formation of a new and small structure be- 
tween the first and second peaks in the RDF. This is asso- 
ciated with 5-fold coordinated atoms and was discussed 
in details in Ref. 30; it is strongest at the lowest tem- 
peratures and essentially disappears at 900 K, consistent 
with the results of Fig. |3] Interestingly, at this tempera- 
ture, the shearing counterbalances the effects of thermal 
disorder, yielding an amorphous state with fewer coor- 
dination defects (lower energy) than the original model. 
The amorphous model is thermodynamically meta-stable 
and the combination of the shearing and thermal effects 
brings the system to a more relaxed structure with fewer 
defects and low energy. 

At 1000 K, now, a completely different behavior is ob- 
served: the RDF is qualitatively different from the low- 
temperature ones, exhibiting much better-defined peaks, 
corresponding in fact to the crystalline state: closer ex- 
amination of the system reveals that it is actually crys- 
tallizing under the yoke of shear. However a complete 
crystallization of the system cannot be reached since it is 
limited by the presence of the walls which remain amor- 
phous during the whole shearing process. 



B. Extent of crystallization 

The extent to which crystallization takes place is best 
understood, and quantified, in terms of an order param- 
eter. For this purpose, we use a topological character- 
ization based on two-ring structures ("blocks") associ- 




FIG. 4: (color online) Fraction of 5-fold coordinated atoms 
at 300 K (a), 600 K (b), 900 K (c) and 1000 K (d). The 
configurations are divided into 12 layers along the y direction; 
layers 1 and 12 are not shown as they correspond to the frozen 
walls. 



ated with local crystalline order f^^EU they are displayed 
in Figj6] (bottom panel). These blocks represent the 
smallest 3D rigid structures that can be extracted from 
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FIG. 5: (color online) Radial distribution function of the sys- 
tem in the non-sheared state and at 12 % strain for the four 
temperatures considered; the curves are shifted for clarity in 
the main diagram and overlaid in the inset. 



the tetrahedral lattice. The first one derives from the 
Wurtzite (hexagonal) structure and contains 12 atoms; 
the second is extracted from the diamond structure and 
is made up of 10 atoms assembled as four 6-fold rings 
back to back. A similar approach was used, for e xamp le, 
in dense-packed systems by Mokshin and Barrat.'^^^ 

We present in Fig. [6] a snapshot of the system under- 
going crystallization; this corresponds to a configuration 
obtained at the end of the simulation at 1000 K and 
a total strain of 12 %. Clearly, the amorphous particles 
(blue) are concentrated near the walls (dark) while atoms 
within a crystalline environment occupy the bulk region 
(yelow). The high concentration of the amorphous par- 
ticles at the interfaces between the walls and the bulk 
region is due to the presence of the walls. Let us men- 
tion that the structure of the walls remain the same dur- 
ing shear deformations. Therefore the crystallization of 
amorphous silicon is limited by the amorphous nature 
of the walls. This situation is similar to that found in 
solid-liquid or amorphous-crystal interface studies. The 
interface is extended over several layers, as for example 
see the Ref. HI 

We present in Fig. [7| the fraction of atoms which can 
be tagged as crystalline within the central region of the 
system as a function of strain for the four temperatures 
considered. At zero strain, about 20 % of all atoms can 
be associated with a wurtzite or diamond structure at 
all temperatures. Since the crystal-like particles are dis- 
tributed uniformly across the simulation box, these are 
associated with local topological fluctuations and do not 
represent a crystalline phase; as discussed in Ref. [4TJ 
the specific fraction of these crystal-like atoms in a well- 
relaxed amorphous network actually depends on the de- 
tails of the interatomic potential. 

Consistent with the results presented in the previous 




FIG. 6: (color online) Snapshot of the system showing the dis- 
tribution of crystalline particles, in light /yellow, at 1000 K; 
amorphous and liquid-like particles are in dark/blue, and the 
rigid walls are in black. The bottom panel shows the blocks 
used to characterize the local crystalline order: Wurtzite 
(hexagonal) structure (right) and diamond structure (left); 
see text for details. 



section, the proportion of crystalline atoms remains ap- 
proximately constant at the three lowest temperatures 
considered (300, 600, 900 K), while fluctuations get larger 
as temperature increases (see inset). These fluctuations 
indicate that, by increasing temperature, a large fraction 
of atoms see their local topological environment change 
between amorphous and crystalline state. The wider dis- 
tribution of the slip planes; slippage is concentrated in a 
few layers at low temperature, as was shown in Fig. [4| 
and involves a larger fraction of the box as the tempera- 
ture rises. 

At 1000 K, however, we observe a qualitative change. 
The fraction of crystalline particles increases steadily — 
slowly for strains < 2 %, then more rapidly up to ~10 %, 
where it reaches a plateau close to 80 %. Such a three- 
stage curve is similar to that observed in the crystal- 
lization of silicon from the liquid phase (see for exam- 
ple Fig. 4 of Ref. SU), suggesting a strain-induced nucle- 
ation process leading to crystallization, limited only by 
the presence of amorphous walls and grain boundaries. 
Complete crystallization however cannot be expected, es- 
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FIG. 7: (color online) Fraction of crystalline particles versus 
strain for the four temperatures considered. 



IV. DISCUSSION 

While several experiments have been conducted to un- 
derstand the micro-mechanisms of plasticity in metallic 
glasses, precise understanding of the atomic processes 
involved in shear propagation is far from being com- 
plete. In particular, several experiments suggest the 
possibility of deformati on- ind uced crystallization, as well 
as nano-void formation.'^U^^'Ell Likewise, recent work on 
single-component Lennard- Jones systems has shown that 
shearing can lead to crystallization from the disordered 
state.-^^-^ However, this has not yet been confirmed for 
more complex, and more stable, disordered systems such 
as amorphous silicon. Indeed, this material is known to 
be very stable at low temperature; it can nevertheless 
crystallize at temperatures well below the melting point 
and numerous studies have focused on crystallization di- 
rectly from the liquid phase^^ or at t he liquid-crystal 
or amorphous-crystal interfaces.^l^^^^ Here, we have 
shown that, when external strain is imposed by a very 
low shear velocity, three different regimes are observed: 
(i) disordered, (ii) annealed and (iii) crystallized. The 
first regime is systematically observed at low tempera- 
tures or at high shear velocity (see Ref. |30| for more de- 
tails). The second and third regimes are clearly observed 
at 900 and 1000 K, respectively. 

As we have discussed previously,^^ when the shear rate 
is faster than thermal relaxation, the steady-state plastic 
regime involves the creation of additional 5-fold, liquid- 
like coordinated atoms that facilitate the constant re- 
organization of the network; the system is thus "more 
disordered" compared to the initial model. This behav- 
ior has also been characterized for various other silicon 
classical potentials by Talati et al7 



The second regime prevails when thermal effects are 
sufficiently large to counterbalance the shear deforma- 
tions. In this case, the applied strain allows the system 
to sample more efficiently the energy landscape, finding 
new basins with low-energy states that can be reached 
due to the available thermal energy. This is what we 
observe at 900 K for all quantities studied: energy per 
atom (Fig.|2|, defect density (Fig. |3|, radial distribution 
function (Fig. [5| and proportion of crystalline particles 
(Fig. [t]). From local to global properties, shearing at 
this high temperature leads to deeper annealing. Inter- 
estingly, even though considerable rearrangements take 
place, the available thermal energy is not sufficient to 
allow nucleation of the crystal and the phenomenon is 
more akin to standard thermal annealing used for pro- 
ducing "good quality" glasses and amorphous materials. 

At 1000 K, however, the combination of temperature 
and shearing leads to crystallization. The evolution of 
the order parameter shown in Fig. [7| follows closely that 
observed for the nucleation from the liquid phased first, 
small nuclei appear with a slow increase of the order 
parameter; this is followed by the rapid growth of the 
largest nuclei until the grains touch and the system be- 
comes poly-crystalline. Remarkably, the time scale over 
which crystallization takes places from the liquid phase 
is more than a 100 times slower than for shearing; this 
is likely due to the fact that diffusion is much slower in 
the amorphous than in the liquid state, and so crystals 
cannot assemble as quickly. 

The validity of our results can be established to 
some degree by noting that they are consistent with 
the corresponding situation in such systems as amor- 
phous alloys based on iron and nickel. Some experi- 
ments have shown, for example, that an external driv- 
ing force can induce crystallization in jammed sys- 
tems tha t can be considered as models for disordered 
materials.l^^EH Yibrations,Ea3Ej shear oscillations ,1^21111 

and steady sheai!^ have been shown to induce crystalliza- 
tion in granular systems of spherical particles, and shear 
oscillations induce crystallization in colloidal glasses.^ 

We note also that crystallization was not observed by 
Talati et alJ in their MD simulations using both the SW 
and the Tersoff potential; this is likely related to the fact 
that they have considered relatively low temperatures. 
We speculate that the details of the crystallization pro- 
cess — temperature range, shear rate, etc. — are poten- 
tial dependent, but the physics is not. 

In contrast to single-component Lennard- Jones, amor- 
phous silicon is a metastable phase that can be an- 
nealed experimentally without crystallization; this al- 
lows us to observe a much slower crystallization pro- 
cess, emphasizing the similarities with classical nu- 
cleation. This suggests that this phenomenon could 
be more much generic that initially thought and 
adds weight to previous observat ion relating shear and 
temperature. B l i" | i^ | i4rii | 4i|42 l44l49]T^-g relation is not per- 
fect, clearly, and further studies are necessary to fully 
understand how one can integrate external deformations 
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into standard nucleation theory. 

V. CONCLUSION 

In summary, we have investigated the behavior of 
amorphous sihcon subject to slow shear deformations 
using classical MD simulations with the EDIP poten- 
tial. While we found a systematic increase of the PED 
and the defects fraction at low temperatures or at high 
shear velocities, the shear deformations, when induced 
by low shear velocities, lead to three different regimes 
that emerge as a function of temperature: increased dis- 
order for temperatures below o:^ 900 K, enhanced relax- 
ation around 900 K, and crystallization at 1000 K. 
The crystallization of amorphous silicon appears to fol- 
low classical nucleation. While the exact parameter space 



for these three regimes to develop certainly depends on 
the choice of interatomic potential, our results, which are 
consistent with a number of other observations, appear 
to be robust and open a new approach for modifying and 
manipulating glasses, and even controlling the formation 
of crystalline nano-structures inside disordered matrices. 
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